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ABSTRACT 

Using the TORUS radiative transfer code we produce synthetic observations of the 
21 cm neutral hydrogen hue from an SPH simulation of a spiral galaxy. The SPH rep- 
resentation of the galaxy is mapped onto an AMR grid, and a ray tracing method is 
used to calculate 21 cm line emission for lines of sight through the AMR grid in differ- 
ent velocity channels and spatial pixels. The result is a synthetic spectral cube which 
can be directly compared to real observations. We compare our synthetic spectral 
cubes to observations of M31 and M33 and find good agreement, whereby increasing 
velocity channels trace the main disc of the galaxy. The synthetic data also show kinks 
in the velocity across the spiral arms, evidence of non-circular velocities. These are 
still present even when we blur our data to a similar resolution as the observations, 
but largely absent in M31 and M33, indicating those galaxies do not contain signifi- 
cant spiral shocks. Thus the detailed velocity structure of our maps better represent 
previous observations of the grand design spiral M81. 

Key words: methods: numerical - radiative transfer - radio lines: galaxies - galaxies: 
individual (M31, M33) 



1 INTRODUCTION 

The last decade has seen huge advances in computational 
modelling of galaxies. Allied with observational develop- 
ments this provides increased opportunities to study and 
understand the interstellar medium (ISM). Comparing mod- 
els and observations provides valuable insight into the dy- 
namics and evolution of the ISM, in particular translating 
the physics incorporated in the simulations into observable 
features. 

High resolution hydrodynamical and magneto- 
hydrodynamical models can now begin to examine how the 
interstellar medium is swept into giant molecular clouds, 
and ultimately forms stars (|Wada Norman|1999 Shetty & 



nearby spiral galaxies M31 (Chemin et al. 2009) and M33 



Ostriker|2006| |Wada|2008l|Dobbs|2008||Tasker Tan|2009| 



Dobbs et al.| 2008| ). The ability of supernovae to trigger 
molecular cloud formation can also be studied using nu- 
merical models. ( |de Avillez Berry||2001| |Dib et"ar]|2006 ). 
At the same time, high resolution surveys now provide 
data on the distribution and properties of molecular clouds 
( Heyer et al.|1998||Engargiola et al.| 2"003| |Rosolows ky|2007 



Brunt et al.| 2009). Advances in observational capability 
also provide increased insight into the characteristics of the 
ISM in external galaxies e.g. recent Hi observations of the 
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( Putman et al.|2009 ) and the irregular galaxies Holmberg II 
(Rhode et al.Ml999) and the LMC (iKim & Park|[2007). 



Additionally, statistical properties of Hi in galaxies can 
be studied using larger samples obtained through surveys 
( [Walter et aTf2008,; ,Tamburro et al.|2009) . Thus comparing 
simulations and observations is particularly timely; however 
in order to provide a direct comparison the output of the 
simulations needs to be processed to produce the same data 
products, such as spectral line datacubes and maps, which 
are obtained from observations. 

An obvious application of such synthetic observations is 
to test whether numerical models can reproduce the struc- 
ture of the interstellar medium seen in observations, thus 
synthetic observations provide an important validation test 
of numerical models. Synthetic observations may also be 
generated with very high spatial resolution compared to ob- 
served data, hence simulations can show what we will view 
in nearby galaxies with future facilities. Moreover, synthetic 
observations are derived from simulated objects with known 
properties (e.g density, temperature, velocity) whereas when 
dealing with real observations these properties are calculated 
from the observations. Within the Milky Way we commonly 
rely on radial velocities to determine distances, however for 
simulations, features can readily be located both in velocity 
and cartesian space. Hence working with synthetic obser- 
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vations makes it easier for features in the observations to 
be related to the physical processes from which they result. 
An important effect in this paper is the impact of velocity 
structure on Hi observations. 

Previous comparisons of velocity structure in simulated 
spiral galaxies with observed data range from simply identi- 
fying clouds and plotting their position in velocity space 
("Dobbs et al.l 12006'; 'Baba et al.||2009|, to sophisticated 



radiative transfer models ( Chakrabarti Whitney^ ^2009^ 
Narayanan et al. 2009) and synthetic galactic plane sur- 
veys (IGomez k CoxpOOi) [Douglas et al.|2010t . |Dobbs et al 



( 2006 ) plotted the distribution of molecular clouds in veloc- 



ity space to show that they reproduce the distribution in 
the Outer Galaxy reasonably well. Baba et al. (2009|) plot- 
ted the location of cold gas using the velocities of the gas in 
their models and an assumed rotation curve, and highlighted 
the disparity with the actual location of the gas, illustrat- 
ing the so called 'finger of God' effect. On larger scales, 



Narayanan et al. (2009 ) and Chakrabarti Whitney (2009) 



perform analysis on SPH (smoothed particle hydrodynam- 
ics) simulations of interacting galaxies, using radiative trans- 
fer codes to determine CO emission and spectral energy dis- 



tributions respectively. Gomez & Cox (2004) and Douglas 



et al. (2010) generate synthetic galactic plane surveys and 



compare their velocity structure to real observations from 
the Leiden/Dwingeloo Hi survey and the Canadian Galac- 
tic Plane Survey respectively. 

On extragalactic scales H i is a powerful tracer used to 
analyse galactic structure, e.g. tracing spiral arms and tidal 
tails, and showing holes and cavities in the ISM. Synthetic 



Hi maps have been constructed previously ( Wada et aL| 
2000||Dib k Burkert|20"05t , though for irregular rather than 
spiral galaxies. These studies compared synthetic Hi maps 
with data from Holmberg H and the LMC, concluding that 
the structure in these galaxies is primarily due to a combi- 
nation of turbulence and thermal and gravitational instabil- 
ities, with super novae contributing to a lesser degree. In this 
paper, we produce synthetic H I maps of a simulated grand 
design spiral galaxy (Dob bs et al.|2008| in which the struc- 
ture is instead dominated by spiral density waves, and the 
gas velocities reveal information about the spiral structure. 
Velocity gradients reveal streaming motions, due to the pres- 
ence of spiral density waves (or at least stellar spiral arms 



which rotate at lower Q(r) than the stars, see Dobbs et al, 



( 2009t ; |Wada Koda| ( |2004| )) which produce spiral shocks 
in the gas. 

As our model galaxy is a perfectly isolated, grand design 
system, which has no effects from interactions with compan- 
ions or high velocity clouds, we expect to see only effects due 



to spiral shocks. Visser (1980) first showed kinks in the ve- 



locities along the spiral arms in M81, due to non-circular 
motions, typically with a magnitude of a few lO's of km s~^ 
(Adler & Westpfahl||1996') . We expect to see similar effects 
in our synthetic observations and our idealised model sys- 
tem allows us to study these effects without the additional 
complexity introduced by environmental interactions. 

We used an SPH code to model the galaxy (unlike the 



grid based calculations by Dib & Burkert (2005) and Wada 
et al.| (|2000)) which we then combine with TORUS, a grid- 



based radiative transfer code. Thus we have the further com- 
plication of converting between SPH and a grid code. We be- 
gin in Section [2] by describing the method used to generate 



synthetic observations, including the SPH to grid conver- 
sion. The technique is used to produce spectral cubes for 
galaxies with orientations like M31 and M33, which are as- 
sessed and compared to real observations of M31 and M33, 
in Section [3] as a means of validating the method. Although 
these galaxies do not display grand design structure to the 
extent of M81 (indeed M33 is a fiocculent spiral), we have 
access to high resolution Hi data for these nearby systems. 



2 METHOD 

This section describes the method used to generate the syn- 
thetic observations, starting from the SPH simulation of the 
spiral galaxy and going through to the production of a 21 cm 
spectral line cube. 



2.1 SPH simulation of the spiral galaxy 

2.1.1 Calculations 

The calculation used in this paper was set up as described 
in Dobbs et al. (2008), but we also provide a description 



below. The calculation constitutes one of a series investigat- 
ing structure and the formation of molecular clouds in spiral 
galaxies, and most notably, includes a full thermodynami- 
cal model of the ISM, and the conversion between atomic 
hydrogen and molecular hydrogen, necessary for producing 
our synthetic observations. We did not however include self 
gravity or magnetic fields, although these have been the sub- 
ject of previous papers ( |Dobbs||2i[)08] [Dobbs k Price||2008 ). 
The presence of magnetic fields is expected to inhibit the 
formation of structure in the galactic disc and reduce the 



strength of spiral shocks (Roberts & Yuan 1970 Dobbs & 
Price|2008| ), resulting in a smoother distribution of gas, and 
spiral arms which are wider and less dense. Conversely, self- 
gravity will increase density in regions which are already 
overdense (enhancing the formation of structure in the gas) 
and increase the likelihood that gas in dense regions will be 
fully molecular. 

The calculation used as the basis for our synthetic ob- 
servations was performed using smoothed particle hydrody- 
namics (SPH), a Lagrangian fiuids method. We modelled 
a gaseous disc, whilst the stellar component was included 
as an external potential. This external potential includes a 
stellar disc and halo, and incorporates a 4 armed spiral com- 
ponent from Cox & Gomez (2002^) with a pattern speed of 

2 X 10 ~^ rad yr~ and a pitch angle of 15"^. The amplitude 

2 



of the stellar spiral perturbation is 1.1 x 10^^ cm^ s 

We did not model the whole disc, but restricted the do- 
main to radii between 5 and 10 kpc, primarily to increase 
resolution. The gas particles were initially distributed ran- 
domly, with a scale height of 400 pc and a temperature of 
7000 K. The particles were assigned circular velocities ac- 
cording to the disc potential, with a 6 km s~^ velocity dis- 
persion. 

The surface density of the galaxy was 10 Mq pc^, in- 
cluding helium, which is comparable to the average surface 
density at the solar radius ( Wolfire et al.|[2003 ). We used 8 
million particles in the calculation thus giving a resolution 
of 500 M0 per particle. 

We followed the thermal evolution of the gas using a 



Synthetic observations of a spiral galaxy 3 



model for the heating and cooUng of atomic and molecular 



Carlo method of Lucy (1999), and can also generate spec- 



gas taken from Glover & Mac Low ( 2007 ) . This model incor- 



porated many processes, including cooling from fine struc- 
ture emission, gas-grain cooling, and heating from photoelec- 
tric emission and cosmic ray ionisation of atomic hydrogen. 

We also included the formation of molecular hydrogen, 
as well as CO. We evolve the abundance of molecular gas 
according to a prescription given by Bergin et al. ('2004 ), and 



provide a full description in D obbs (2008 ) . We assign each 
particle a molecular gas fraction, and calculate the rate of 
formation of H2 on grains, and the rate of destruction of H2 
by photodissociation and cosmic ray ionization, to update 
the molecular gas fraction after each time step. We noted in 
our previous paper that determining the photodissociation 
rate requires an estimate of the column density in order 
to take into account self shielding. As described in Dobbs" 
( [2OO 8 ) we adopt a constant length of 35 pc to calculate the 
column density, which is the average distance to a BO star 
in the Milky Way. Given the H2 fraction for a particular 
particle, we determine the Hi density used for this paper 
from n{H) = n — 2n(if2), where n{H) , n{H2) and n are the 
number densities of Hi, H2 and the total number density. 



2.1.2 Evolution of galaxy disc 



The evolution of the gas disc is described fully in |Dobbs et al.| 
( 2008 ) but we also provide a brief summary here. As the disc 



evolves, the spiral pattern emerges and narrow spiral arms 
develop. Most of the gas cools from the initial temperature 
of 7000 K and typically 70 % of the gas (and all the gas in 
the spiral arms) is < 150 K. Clumps of cold gas accumu- 
late into more massive clouds as they pass through a spiral 
shock. The molecular gas (which constitutes around 30 % 
of the total gas) lies predominantly in the spiral arms, and 
typically forms from the atomic medium over time scales of 
order 10 Myr. As dense clumps leave the spiral arms, they 
are sheared into short spurs. Some cold clumps, or spurs, 
are able to survive between the arms, though they do not 
tend to contain much molecular gas. We do not have ob- 
servational measurements of the proportions of warm and 
cold gas in external galaxies. The fractions in our simula- 
tions are not dissimilar from the solar neighbourhood, the 
only region for which there is an observational indication of 
the fractions of cold and warm HI dHeiles k Troland|2003| ). 
In the absence of observational evidence to the contrary we 
can only assume that other nearby large spiral galaxies are 
not vastly different in terms of the amount of gas in different 
phases. 

We run the simulation for a total of 320 Myr. For the 
comparisons in this paper, we take the output at a time of 
250 Myr, by which point the distribution of gas in differ- 
ent phases, and the amount of molecular gas, has reached a 
roughly steady state. 



2.2 AMR grid construction 

We use the radiative transfer code TORUS ( Harries||2000 ) to 
generate synthetic Hi observations. TORUS is a grid-based 
radiative transfer code that uses Adaptive Mesh Refine- 
ment (AMR) to provide variable spatial resolution. TORUS 
can perform radiative transfer calculations, using the Monte 



tral energy distributions, images and spectral cubes. The 
code has frequently been applied to models of stellar discs 
and performs well in benchmark tests, even at high optical 
depths ( |Pinte et al.]|2009| ). 

To generate a synthetic Hi data cube we first need to 
convert from the particle representation of SPH to the AMR 
grid representation of TORUS. The TORUS AMR grid is con- 
structed using the octree method in which the grid initially 
comprises eight cells (one octal) with two cells in each spa- 
tial dimension. The grid is refined by specifying a condi- 
tion which determines when a cell is split into a further 
eight cells. For this calculation the grid cells were split if 
the mass of the cell exceeded a given limit (mass per cell 
limit). A mass per cell limit gives higher spatial resolution 
in regions of high density, which is similar to the effective 
spatial resolution of the SPH method. A maximum Hi mass 
per ceU of 2.5 x 10^^ g (1260 M©) was used to split the grid, 
which resulted in an AMR grid comprising 678815 octals and 
4751706 unique cells. This mass per cell limit corresponds 
to approximately 9 SPH particles per grid cell (for particles 
with an average atomic hydrogen fraction) and was chosen 
to give sufficient spatial resolution to represent the impor- 
tant features of the model galaxy, within the constraints of 
the available computer memory. Should a calculation require 
higher spatial resolution, the mass per cell limit can be de- 



creased to give smaller cells in the AMR grid ( ^Douglas et al. 



2010) 



Once the grid structure has been generated H i density, 
temperature and velocity values are assigned to each cell 
using values from the SPH particles, smoothed with an ex- 
ponential kernel. Construction of an AMR grid from SPH 
particles, using this method, was tested by |Acreman et ah] 
( |2010J using an azimuthally symmetric circumstellar disc 
benchmark. In the case of a spiral galaxy it is also neces- 
sary to represent accurately structure, such as spiral arms, 
clouds and spurs, which are not present in an azimuthally 
symmetric disc. We found that the method of |Acreman et al^ 
( |2010J gave a good representation of structure within the 
disc and a good representation of the total Hi mass pro- 
vided one modification was made. The method normalises 
density values by the sum of the SPH kernel weights, if the 
sum of the weights exceeds a specified threshold, in order 
to reduce noise in the interior of the distribution. 



Acre manI 

2010) normalise when the sum of the weights exceeds 



et al. 



0.3, but we found this value gave a total Hi mass which 
was too large, by 11 per cent, relative to the Hi mass of the 
SPH particles. The SPH to grid conversion was performed 
with a range of normalisation thresholds, in order to test 
the impact on the representation of the total H i mass, and 
the results are plotted in Fig. ^ A normalisation threshold 
of 0.5 gives a mass which is too small by only 0.3 per cent, 
so we use this value as the threshold above which the SPH 
kernel smooth is normalised. 

The Hi density in the galaxy midplane, as represented 



on the AMR grid, is shown in Fig. 2(a) For purposes of 
comparison a similar plot from the original SPH particles is 
shown in Fig. 2(b) where the particles have been rendered 
using the SPLASH plotting software ( |Price|2007| ). This fig- 
ure shows that the spiral arms are well represented in both 
cases and the SPH particle to AMR grid conversion has been 
effective in preserving these features. Figure [3] again shows 
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Figure 2. Hi density (in g/cm ) in the galaxy midplane, as represented on the TORUS AMR grid (left) and as represented with SPH 
particles (right). The spiral arms are well represented in both cases demonstrating that the conversion from SPH particles to AMR grid 
preserves spiral arm structure. 
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Figure 3. Hi density (in g/cw?) in the galaxy midplane, as represented on the TORUS AMR grid (left) and as represented with SPH 
particles (right). The AMR representation has the grid overplotted. This figure is similar to Fig. |2] but shows a smaller domain to 
emphasize that the smaller scale structures are well represented on the grid. 



the midplane density in the grid and SPH representations, 
but for a reduced domain, with the AMR mesh over plotted 



in Fig. 3(a) to emphasize how the enhanced AMR resolu- 
tion in regions of high density enables good resolution of 
small scale features. 



2.3 Opacity and emissivity calculation 



The relatively simple nature of the Hi 21cm line transi- 
tion means that once the density and temperature of a cell 
are known the emissivity and opacity can be calculated and 
stored on the AMR grid. The opacity is given by 
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Figure 1. Fractional error in Hi mass on the AMR grid, relative 
to the H I mass of the SPH particles, as a function of threshold 
for normalisation of the SPH kernel smoothing. Positive values of 
mass error indicate that the grid mass is too large. The lines are 
line sections joining points and not a best fit line. 



Sc^hAo n(H) 
327Tkiyo Ts 



(1) 



(Rohlfs & Wilson 2004) where c is the speed of light, h is 



Planck's constant, Aq is the Einstein probability emission 
coefficient, k is the Stefan-Boltzmann constant, and is the 
frequency of the Hi transition. The number density of Hi 
(n(H)) and the excitation temperature (Ts) are taken from 
the density and temperature values on the AMR grid. The 
emissivity is calculated from the opacity using Kirchoff's 
law {tu = nuBu where Bu is the Planck formula) and the 
Ray leigh- Jeans approximation giving 



= -^g^ n(H)0(^;) . 



(2) 



The profile function cj) {v) is assumed to be a Gaussian 
with a thermal line width ath given by 



2 

(^th — 



kT 
ruH 



(3) 



where T is the temperature of the grid cell and ttih is the 
mass of a hydrogen atom. Unlike Wada et al. ( 2000| ) we do 
not add a term to represent microturbulence as the major- 
ity of the gas in our simulation is at temperatures where 
we expect the thermal line width to dominate (see fig. 4 of 
Dobbs et al. (2008)). As the contribution of the line proffie 



to the opacity and emissivity depends on the velocity chan- 
nel under consideration the line profile is calculated during 
the radiative transfer step and is not stored on the AMR 
grid. 



2.4 Radiative transfer calculation 

Once the AMR grid is set up a spectral cube is generated by 
carrying out ray tracing operations through the AMR grid. 



The method is similar to that of Rundle et al. ( 2010[ ) and is 
summarised below. 

A ray is traced through the grid, starting from the ob- 
server's position and proceeding along the ray direction to 



distances further from the observer. The contribution of the 
grid cell to the pixel brightness is calculated using the emis- 
sivity and opacity taken from the AMR grid. Each cell con- 
tributes to the intensity in the pixel but with absorption by 
the accumulated optical depth between that cell and the ob- 
server. There is an iterative procedure in which the intensity 
of the pixel is updated from an old value to a new value 
II with the contribution of a grid cell according to 



= + (l 



(4) 



where Su = ^ ^ dr \s the optical depth of the cell and rtotai 
is the total optical depth between the cell and the observer. 

The model grid is rotated to represent a galaxy with 
the required position angle and inclination angle. A number 
of parallel rays, one per pixel, are traced through the grid 
to calculate emission for each pixel in each velocity channel 
of the cube. The ray trace through each individual cell is 
sub-divided into smaller steps to give better sampling of the 
velocity gradient across the cell. If the velocity difference 
across the cell, projected onto the direction of the ray, is 
then the ray trace is decomposed into 5^ separate steps. 

The synthetic observations are converted to units of 
brightness temperature, using the Ray leigh- Jeans approx- 

■V 2 

imation, by multiplying by where A is the wavelength 
of the transition. Continuum emission from the cosmic mi- 
crowave background is included as a blackbody component 
with a temperature of 2.73 K. 



3 RESULTS 

Our technique was used to generate synthetic spectral line 
data cubes resembling radio observations of M31 and M33, 
two Local Group spiral galaxies, for which arcminute-scale 
21 cm observations were available. The model galaxy was 
not set up to model M31 or M33 specifically, so although 
we expect the synthetic data to reproduce general features 
seen in Hi observations of external galaxies we do not make 
a detailed quantitative comparison with the observations. 
However comparisons of the synthetic data with observed 
data do provide a valuable validation test of the method. 



3.1 M31 

The model galaxy was given an inclination angle of 77.5 de- 
grees and a major axis position angle of 220.0 degrees (mea- 
sured clockwise from North) in order to match the observed 
orientation of M31 in equatorial co-ordinates (see figure 9 
(|2009)). The position angle of the model 



of Chemin et al, 



galaxy is offset by 180 degrees so that it rotates in the same 
direction as M31. A velocity offset of —300 km/s was ap- 
plied to the data cube velocity channels consistent with the 
systemic velocity determined by Chemin et al.| ( |2009 ). The 
orientation of the galaxy can be seen in the column density 
plot in Fig. 4(a)[ The spectral cube was constructed using 
600 spatial pixels of size 10^° cm and 250 velocity channels 
over a range of 840 km/s. The velocity channels were chosen 
to give comparable velocity resolution to the observations of 
Chemin et all (120091). 



The emission from the synthetic data cubes is plotted as 
"renzograms" (Fig. 4(b) and Fig. 4(c) ) which plot a single 
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(c) Simulated data blurred with 3 pixel Gaussian (d) Observed data 

Figure 4. Top left: column density from the model galaxy orientated to be like M31. Top right: contours of constant brightness 
temperature in different velocity channels (renzogram) overlaid on summed intensity for the model galaxy. There is one contour per 
velocity channel with different velocity channels indicated by different colours (—500 km/s for the red contour and —100 km/s for the 
black contour). The data have been blurred with a one pixel Gaussian. Lower left: as for top right but blurred with a three pixel Gaussian. 
Lower right: renzogram plot from the M31 observation of [Chemin et al.|2009| 
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Figure 5. A spatially averaged line profile from the simulated 
M31 galaxy (solid line, left side y-axis). The average is taken 
over the whole data cube and background subtracted. A similar 
spatially averaged line profile from the observations of [Uhemin| 
|et al. 120091 is ^^^^ plotted (dashed line, right side y-axis). 



Figure 6. Line profiles at a point from M31 synthetic data with 
1000 velocity channels (solid line and squares) and 250 velocity 
channels (dashed line and crosses.) 



level contour map for several velocity channels. The con- 
tours are of constant brightness temperature with different 
velocity channels shown in different colours. The contours 
are at 50 km/s intervals starting at —500 km/s (red con- 
tour) and ending at —100 km/s (black contour). The con- 
tours are overlaid on a grey scale plot of brightness tempera- 
ture summed over all velocity channels. These synthetic data 
cubes have been spatially blurred with a Gaussian filter of 
1 pixel (Fig. |4(b")] ) and 3 pixels (Fig. 4(c) ) width to simulate 
the effect of observing with an instrument with finite spatial 
resolution (1 pixel corresponds to 0.14 arcmin and 3 pixels 
corresponds to 0.43 arcmin at a distance of 785 kpc). A simi- 
lar plot of the observational data presented by |Chemin et al.| 
(2009) is shown in Fig. 4(d) As the rotation velocity of the 
model galaxy (220 km s~ ) and the rotation velocity of M31 
(230 km/s in the outer galaxy Chemin et al.^2009j are sim- 



ilar, the velocity channels which are contoured in Fig. 4(d) 
are the same as in Fig. 4(b) and Fig. |4(c)| 



A good general agreement is seen between the synthetic 
and real data whereby the galaxy's main outer disc is traced 
as the velocity channels increase. The synthetic data with 
the least blurring (i.e. the highest spatial resolution) show 
structure associated with the spiral shock in the disc. This 
structure is not particularly evident in the observed data 
for M31. This may be because of lack of resolution (it is 
also more difficult to distinguish the spiral structure in our 
3 pixel blurred image) or simply because the spiral shocks 
are too weak or nonexistent. In fact Efremov (2009) argues 



that there is a spiral shock in one arm, but not the other, 
based on observations of stellar gradients across spiral arms 
and the distribution of star complexes, so the presence of 
spiral shocks in M31 is rather ambiguous. 

A line profile from the synthetic data, blurred with a one 
pixel Gaussian and spatially summed over the whole galaxy, 
is plotted in Fig. [s] (solid line). A similar observed global line 
profile from the data of Chemin et al. ( 2009 ) is also plotted 



profile can be quantitatively compared to the observed line 
profile. Both profiles show peaks at around —100 km s~^ and 
—500 km s~^ although M31 has a slightly greater circular 
velocity, so the peaks are slightly further apart. The ob- 
served profile is more extended than the synthetic profile, in 
that the emission falls away less rapidly towards -600 km/s 
and zero; the effect is more pronounced towards more nega- 
tive velocities. Figure 7 of | Chemin et al.] |2009| shows that 
material with velocities around -600 km/s is to be found 
at approximately one third of the disc outer radius (0.8 de- 
grees from the centre of M31). This material would lie within 
the 5 kpc inner radius of our model galaxy and would not 
be represented. The observed profile also shows asymmet- 
ric structure between the peaks, which is not present in the 
synthetic line profile, which we also expect to come from the 
central regions of the galaxy. There is likely to be complex 
structure present in the inner regions of a real galaxy, which 
is not represented in our model galaxy, and this structure 
can noticeably affect the global line profile. 

Many of the line profiles from individual pixels contain 
structure which is better represented in data with higher 
velocity resolution. The synthetic data cube was regenerated 
using 1000 velocity bins, over the same velocity range, and 
the effect on an example line profile can be seen in Fig. |6] 
(solid line with squares is 1000 velocity channels, dashed 
line with crosses is 250 velocity channels). This profile is 
taken from x — —1.59 kpc, y — 3.91 kpc which is in one 
of the spiral arms. The line profile from the data cube with 
250 channels shows that the shape of the line is much less 
accurately represented than with 1000 channels and some 
features of the line (e.g between —165 and —170 km/s are 
not represented in the lower resolution case. The presence 
of structure with intrinsically narrow velocity widths was 



(dashed line). As the circular velocity of the model galaxy 
is similar to the rotation velocity of M31 the modelled line 



observed in VLA observations of M31 by Braun (1990) and 
our model galaxy exhibits similar narrow velocity widths in 
its spectral features. For accurate quantitative work the data 
which are output from the radiative transfer code must be 
represented with sufficient velocity resolution to represent 
accurately these narrow features of the line profiles. 
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3.2 M33 



The model galaxy was given an inclination angle of 50.0 de- 
grees and a major axis position angle of 20.0 degrees (mea- 
sured clockwise from North) in order to match the observed 
orientation of M33 in equatorial co-ordinates. These values 
are representative of the Hi distribution fitted by |Corbelli1 
& Schneider (1997). A velocity offset of —180 km/s was ap- 
plied to the data cube velocity channels in order to match 
the observed approach velocity( Corbelli Schneider|1997 ). 
The cube comprises 600 spatial pixels of size 10 cm and 
500 velocity channels over a range of 360 km/s. 
A column density plot is shown in Fig. ' 



7(a 



which shows 

a more face-on view than M31 (see Fig. 4(a)). Renzogram 
plots, overlaid on a grey scale plot of summed intensity, are 



shown in Fig. 7(b) and Fig. 7(c) for synthetic data blurred 
with a 1 pixel and 6 pixel Gaussian respectively (1 pixel 
corresponds to 0.15 arcmin and 6 pixels corresponds to 0.92 
arcmin at a distance of 730 kpc). The minimum velocity is 
—342 km/s and contours increase in steps of 36 km/s to a 
maximum of —18 km/s. A corresponding contour plot from 
the observational data of Putman et al. (2009) is shown 



in Fig. 7(d) As M33 has a smaller rotation velocity than 
our model galaxy (100 km/s compared to 220 km/s) the 
velocities which are contoured in Fig 7(d) have been scaled 
by a factor of 100/220 compared to the velocities contoured 
in the synthetic data plots. For both the synthetic data and 
the observed data the contour interval is 16 per cent of the 
rotation velocity. 

We again get a good qualitative comparison between the 
synthetic and real observations, taking into account that we 
do not model the inner 5 kpc of the galaxy. The general be- 
haviour of the contours is as expected i.e. the main disc of 
the galaxy is traced out as velocity increases and there are 
perturbations in the velocity contours due to non-circular 



velocities. Figure 7(b) clearly shows perturbations in the 
contours associated with the spiral arms whereas lower res- 



olution data (Fig. 7(c) ) show much less detail but exhibit a 
broadening of the contours in the region of the spiral arms. 
Our simulated galaxy is a grand design spiral which shows 
non-circular motions due to the presence of spiral shocks. In 
contrast M33 is a flocculent spiral with a relatively weak spi- 
ral perturbation. The structure in M33 is instead largely due 
to gravitational instabilities and supernovae feedback which 
are not included in the model galaxy. Moreover even in the 
absence of spiral shocks there are likely to be non-circular 
motions due to other mechanisms (e.g. tidal interactions and 
warps). The simulated data show what would be observed 
in an idealised case where a strong spiral perturbation dom- 
inates the structure of the galaxy. 

3.3 Comparison to optically thin limit 

In the optically thin limit the absorption term in the radia- 
tive transfer equation can be neglected i.e. 

dl 



ds 



(5) 



and the intensity in an image pixel can be found by simply 
integrating the emissivity along the required line of sight. In 
this case the intensity is proportional to the column density. 
As a validation test of the method the ratio of intensity from 
the radiative transfer code to the intensity derived from the 




-300 -290 -280 

Velocity channel (km/s) 



Figure 9. Line profiles from x = —0.972 kpc, y = —0.3626 kpc. 
The solid line is from a calculation using 250 velocity bins, the 
dashed line is from a calculation using 1000 velocity bins. The 
velocity integrated emission is overestimated in the 250 bin profile 
due to a lack of resolution in velocity space. 



optically thin approximation was calculated for each pixel 
in the data cube. The ratio of these intensities should be 
unity when the optical depth is small, and the optically thin 
approximation holds, and should decrease below unity as 
the optical depth increases. 

Figure [S] plots the ratio of intensity from the radiative 
transfer calculation to intensity from the optically thin ap- 
proximation against optical depth for the simulated M31. 



Figure 8(a) is from the results presented in the previous sec- 
tion with 250 velocity bins over a range of 840 km/s. Pixels 
in which the column density is less than 10^° cm~^ are not 
plotted in order to exclude pixels which are not associated 
with the galaxy. The intensity ratio at low optical depths is 
close to unity and there is a trend of decreasing intensity, rel- 
ative to the optically thin limit, as optical depth increases, as 
expected. However there is a considerable amount of scatter 
and there are some points with a ratio greater than unity, 
up to a maximum value of 1.21. An example line profile from 
a pixel with a ratio greater than unity is shown in Fig. [9] 
This pixel is located at x = —0.972 kpc, y = —0.3626 kpc 
and is not included in Fig. [S] as the column density is be- 
low the cutoff threshold. However it is used here as a clear 
example of how the velocity integrated emission in a pixel 
can be overestimated due to a lack of resolution in velocity 
space. Figure 8(b) shows the same plot as Fig, 8(a) but for 



the synthetic data cube with 1000 velocity bins. The points 
with ratios greater than unity are absent and there is less 
scatter on the distribution. The velocity resolution clearly 
has a significant impact on the results of this validation test. 

A similar plot showing the ratio of the intensity from the 
radiative transfer calculation to intensity from the optically 
thin approximation against optical depth for the simulated 
M33 is shown in Fig. |10| As before pixels with a column 
density less than 10^*^ cm~^ are not plotted. The velocity 
resolution of the synthetic M33 cube is high enough that 
the comparison with the optically thin limit is good and 
does not show the scatter seen in the M31 cube when only 
250 velocity channels are output. 
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X (kpc) 



,h3gm 35m 33m 32m 

Right Ascension (J2000) 

(c) Simulated data blurred with 6 pixel Gaussian (d) Observed data 

Figure 7. Top left: column density from the model galaxy orientated to be like M33. Top right: contours of constant brightness 
temperature in different velocity channels (renzogram) overlaid on summed intensity for the model galaxy. There is one contour per 
velocity channel with different velocity channels indicated by different colours. The data have been blurred with a one pixel Gaussian. 
Lower left: as for top right but blurred with a six pixel Gaussian. Lower right: renzogram plot from the M33 observation of Put man et alT] 
[20091 



3.4 Velocity perturbations in renzograms 

We develop a simple analytical representation of the ve- 
locity perturbations seen in the renzograms, in order to 
gain a quantitative understanding of how changes in veloc- 
ity within the galaxy relate to perturbations in renzogram 
contours. The following analysis uses a cylindrical polar co- 



ordinate system with its origin at the centre of the galaxy. 
The r co-ordinate represents the distance from the galactic 
centre in the midplane, the z co-ordinate is the distance out 
of the midplane and 6 is the cylindrical polar angle. It is 
assumed that all gas moves in the azimuthal direction only. 

The line of sight velocity, towards the observer, of a 
parcel of gas is 
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(a) (b) 

Figure 8. Ratio of intensity from radiative transfer calculation to intensity from the optically thin approximation, plotted against optical 
depth (for the simulated M31). Left: from a calculation using 250 velocity bins. Right from a calculation using 1000 velocity bins. 



o 




Figure 10. Ratio of intensity from radiative transfer calculation 
to intensity from the optically thin approximation plotted against 
optical depth for the simulated M33. 



Vt COS sin i + v^ys 



(6) 



where vt is the magnitude of the tangential velocity (velocity 
in the azimuthal direction) for a given parcel of gas, i is the 
inclination of the galaxy and ^;sys is the systemic velocity 
of the galaxy. For an observation with the systemic velocity 
subtracted this gas will be seen in velocity channel 



'^^ch = '^^t cos sin i (7) 

In the absence of any deviation from a constant circular 
velocity the renzogram contours have a constant value of 0. 
However, the effect of the spiral shock means that in a given 
velocity channel we expect to see both gas moving at the 
circular velocity of the galaxy {vt = Vc) and other gas which 
has just passed through a spiral shock (post-shock gas). This 
causes an azimuthal perturbation in the renzogram contour. 

We assume that the effect of passing through the shock 
is to reduce the speed of this gas to vt = fvc (0 < / < 1) 
without changing its direction. Gas moving at the circular 
velocity will be seen at a position 

COSt/c = 

Vc sm i 

whereas post-shock gas is seen at a position 

cos Ops = -r^. — : (9) 

jVc smz 

The linear separation in the plane of the disc {d) of circular 
and post-shock gas is given by 



(8) 



2r sin 



f Oc — Ops \ 

\ 2 j 



(10) 



which can be re-written, using trigonometric identities and 
substituting from eqn. [8] and eqn. [5] as 



where 



(1 



2x1/2 X 




(11) 



X = (12) 

Vc sm I 

When viewed on the plane of the sky this distance is reduced, 
due to projection effects, by a factor of 
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(sin^ + cos^ i cos^ 0^ "^^^ = (l + (cos^ i — l) x 



2\l/2 



(13) 



Hence the angular size of the perturbation in the velocity 
contour (Aa) is 



Aa 



2r 



gal 



X (l + (c( 




,1/2 



(14) 



where r is the distance of the perturbation from the galactic 
centre and i^gai is the distance to the galaxy. 

At X = (i.e. the channel corresponding to the sys- 
temic velocity) there is no velocity perturbation according 
to eqn. |14| In this channel gas which is moving in circular 
motion has no line of sight component. As we have assumed 
that the velocity perturbation is simply a reduction of the 
circular motion there is no mechanism for changing the line 
of sight component. In the renzograms derived from the syn- 
thetic observations there are perturbations seen in this ve- 
locity channel indicating that the velocity perturbation in 
our model galaxy is not simply a reduction of the circular 
motion but has another component which affects the line of 
sight velocity. 

The angular size of the velocity perturbations seen in 

against 

X in Fig. 



the renzograms (as predicted by eqn 14) is plotted 

for four inclination angles; degrees (solid 
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line), 30 degrees (dashed line), 60 degrees (dot-dashed line) 
and 90 degrees (dotted line). This plot assumes the galaxy 
is at a distance of 785 kpc (i.e the distance to M31) with a 
shock 7 kpc from the galactic centre (representative of our 
model galaxy). The shock is assumed to reduce the circular 
velocity by a factor / = 0.9 (i.e. a perturbation of 20 km/s 
on a 200 km/s circular velocity). A similar plot showing the 



size of the velocity features in kpc is shown in Fig. 11(b) 



The predicted size of the features is in agreement with the 
features seen in the renzograms of the synthetic observa- 
tions (i.e. a typical size of approximately a kpc). In all cases 
the size of the velocity perturbation tends towards zero as x 
tends towards zero, due to the assumption that the pertur- 
bation only affects the circular velocity (as discussed above). 
As X increases towards / the size of the feature increases be- 
cause the circular velocity has a larger component in the line 
of sight. The renzograms for our synthetic observations also 
show larger velocity kinks in channels further from the sys- 
temic velocity. There is a change in the sign of Aa when x 
becomes negative (eqn 14 is an odd function of x) due to the 
assumption that the shock always reduces the circular ve- 
locity. This results in post-shock gas always moving towards 
the velocity channel of the systemic velocity (or zero if the 
systemic velocity has been subtracted). The size of the ve- 
locity kinks increases as the inclination angle increases, due 
to projection effects, however in practice it will be difficult 
to see these features for nearly face-on cases because very 
high velocity resolution would be required (for nearly face-on 
cases the range of x corresponds to a small velocity range). 



4 CONCLUSIONS AND FUTURE WORK 

We have presented a method for generating synthetic spec- 
tral cubes of the 21 cm hydrogen line from SPH simulations 



of spiral galaxies. The method successfully maps the den- 
sity, temperature and velocity from the SPH particles onto 
an AMR grid, while preserving important structures (e.g. 
spiral arms and spurs) and accurately representing the total 
mass. Synthetic data cubes are generated using a ray trac- 
ing method. The synthetic data show good agreement with 
observations of M31 and M33 whereby increasing velocity 
channels trace out the main disc of the galaxy. Velocity con- 
tours of the synthetic data show perturbations due to non- 
circular motions, similar to those already observed in M81 



( |Adler Westpfahl|1996l|Visser|1980l|Rots Shane|1975| ), 

which are not seen in the observations of M31 and M33. 

Our model galaxy is a grand design spiral galaxy and 
shows velocity structure associated with the spiral perturba- 
tion. The method of generating synthetic observations can 
also be applied to simulations in which velocity structure 
is generated by other mechanisms. Internal mechanisms, 
such as self- gravity and stellar feedback, can generate veloc- 
ity structure, and indeed are dominant in flocculent spiral 
galaxies. External influences can also affect the Hi morphol- 
ogy of a spiral galaxy e.g. interactions with a companion 
galaxy, high velocity clouds or the intracluster medium. As 
the majority of galaxies reside in groups or clusters, it is 
likely that external environmental effects will influence the 
Hi structure of most real galaxies; our model galaxy is per- 
fectly isolated and our results show idealised behaviour in 
the absence of external influences. 

A large parameter space could potentially be studied, 
involving galaxies with different internally and externally 
generated velocity structure. Some of these aspects can al- 
ready be modelled (e.g. a galaxy undergoing an interaction. 



|Dobbs et al.|(2009| ) and some require further model devel- 
opment (e.g inclusion of stellar feedback, which is ongoing). 
A large number of models must be be run in order to gen- 
erate the SPH input for the radiative transfer calculation. 
Once such a library of simulations was available, one could 
compare these with real observations in order to understand 
the velocity structure of spiral galaxies. However given the 
potentially complex nature of the velocity structures gen- 
erated there would need to be a reliable way of matching 
an observed galaxy with its counterpart in the synthetic 
observation library. In our present models, we can relate 
the strength of the spiral shock, and the inclination of the 
galaxy, to the size of the perturbations in the renzograms. 
When we include additional physics, e.g. self gravity and 
feedback, we can investigate whether these perturbations 
are still observable (for a given shock strength) or whether 
they are overwhelmed by other motions in the gas. 

Our method has also recently been applied to an ob- 
server placed inside the galaxy ( Douglas et al.|20To l) to gen- 
erate a synthetic galactic plane survey. Hi self absorption 
features and the conversion of cold H i to molecular clouds. 



as seen in a real galactic plane survey (e.g. Taylor et al. 
( |2003| ); [StIi et al.| ( [2006t ), are seen in the synthetic data. A 
powerful future application of using a simulation for gener- 
ating a survey is that given a series of time frames we will be 
able to trace the evolution of molecular clouds, and related 
observed features, to the physical properties of the material 
in the cloud. The method can also be extended to generate 
synthetic observations of other tracer species, such as CO, or 
dust, for simulated surveys or simulated external galaxies. 



12 David M. Acreman, Kevin A. Douglas, Clare L. Dohhs, Christopher M. Brunt 






(a) (b) 

Figure 11. Left: Angular size (in arcmin) of velocity perturbations, according to eqn. |14[ for a galaxy at 785 kpc with a shock at 7 kpc 
from the galactic centre. The shock strength parameter is / = 0.9. Four inclinations are shown: degrees (solid line), 30 degrees (dashed 
line), 60 degrees (dot-dashed line) and 90 degrees (dotted line). Right: as for previous figure but with sizes in kpc. The parameter x (see 
eqn|12|) is the velocity channel scaled by the circular velocity projected onto the line of sight (-Ucsini). 
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